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Viscosity in cosmological simulations of clusters of galaxies 
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ABSTRACT 

The physics of the intracluster medium, in particular the values for the ther¬ 
mal conductivity and the viscosity are largely unknown and subject to an ongoing 
debate. Here, we study the effect of viscosity on the thermal state of the intr¬ 
acluster medium using three-dimensional cosmological simulations of structure 
formation. It is shown that viscosity, provided it is not too far off from the 
unmagnetised Spitzer value, has a significant effect on cluster profiles. In par¬ 
ticular, it aids in heating the cool cores of clusters. The central cooling time of 
the most massive clusters in our simulation is increased by more than an order 
of magnitude. In large clusters, viscous heating may help to establish an entropy 
floor and to prevent a cooling catastrophe. 

Subject headings: galaxies: active - galaxies: clusters: cooling flows - X-rays: 
galaxies 


1. Introduction 

Cooling by bremsstrahlung and line emission leads to a loss of pressure support in the 
centers of galaxy clusters, which, in the absence of non-gravitational heating, would cause a 
slow, subsonic inflow of gas towards the center of the gravitational well. As a result, cluster 
cores should cool and accrete gas at rates of hundreds and more solar masses per year. This 
scenario is in conflict with observational evidence that indicates that mass deposition rates 
are consistently lower than predicted. Moreover, the gas temperatures in cluster centers are 
maintained typically above ~ 2 keV (Peterson et ah 2001). 


More evidence for non-gravitational heating in clusters comes from cluster scaling rela¬ 
tions. These relations show departures from self-similarity: In the absence of non-gravitational 
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heating and radiative cooling, the entropy is expected to scale linearly with the mean cluster 
temperature. However, observations by Ponman, Sanderson & Finoguenov (2003), Pratt & 
Arnaud (2005) and Piffaretti et al. (2005) indicate a scaling of entropy roughly according 
to Moreover, they reveal a systematic excess of entropy in low-mass clusters (e.g., 

Ponman, Sanderson & Finoguenov 2003). Churazov et al. (2001); Briiggen et al. (2002); 
Bruggen & Kaiser (2002) have argued that heating by a central AGN can keep the ICM 
from cooling dramatically in the center. 


Thermal conduction has been put forward to explain the absence of soft X-rays from 
galaxy clusters (Kim & Narayan (2003) and references therein). Cosmological simulations 
with thermal conduction have been performed by Dolag et al. (2004) and Jubelgas et al. 
(2004). If thermal conduction is at work, other transport processes such as viscosity are 
bound to be important, too. Thermal conduction transports energy and is mediated mainly 
by the faster electrons. Viscosity, on the other hand, transports momentum and is mediated 
primarily by the more massive ions. It is not clear that the suppression factors of both 
transport processes should be the same. Whether they are the same may depend on the 
scale magnetic fluctuations extend to. This scale may be much larger than the gyroradii of 
electrons and ions (in which case suppression factors could be comparable) or it could be 
comparable to the ion gyroradius. The magnitude of the suppression factor is motivated by 
various theoretical arguments (e.g., given in Narayan & Medvedev 2001). However, we note 
that the precise value of the suppression factor is highly uncertain and, depending on the 
nature of magnetic turbulence, may even exceed the Spitzer value (Cho et al. 2003) or be 
supressed well below it. 


Based on observations of the Perseus cluster, it has been suggested by Fabian et al. 
(2003) that viscosity may play an important role in dissipating energy injected by the central 
AGN. The case for this is based on the existence of long, straight Ha-hlaments that appear 
to rule out the presence of strong turbulence in the cores of galaxy clusters. The Reynolds 
number for a fluid flow whose viscosity is suppressed with a factor, /, with respect to the 
Spitzer value, is given by 


Re ~ 50 /o .3 ne-3 Tkpc usoo , (1) 

where is the electron number density in units of 10“^ cm“^, Lkpc the typical size of an 
eddy in units of kpc, Usoo the associated velocity in units of 500 km s“^ and T7 the temper¬ 
ature of the fluid in 10^ K. As such values of Reynolds numbers are below the critical value 
that separates laminar and turbulent regimes, it was concluded that viscosity can play an 
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appreciable role in the ICM, provided that viscosity is not heavily suppressed. 


Subsequently, heating by viscous dissipation of AGN-induced motions has been simu¬ 
lated by Ruszkowski et ah (2004a,b); Briiggen et ah (2005). It was concluded that, provided 
viscosity is not suppressed signihcantly with respect to its unmagnetised value, viscous dis¬ 
sipation of AGN-induced motions can balance the radiative losses in the IGM. In Reynolds 
et ah (2005) the effect of viscosity on the evolution of radio bubbles was studied, and it 
was found that viscosity had a stabilising effect on underdense bubbles. Fujita et ah (2004) 
have studied the dissipation of motions induced by acoustic-gravity waves in cluster cores. 
They hud that, provided the wave amplitude is large enough, they can suppress the radiative 
cooling of the cores. Kim et ah (2005) have investigated the heating of clusters by dynamical 
friction. They concluded that friction can be an important supplier of heat but is unlikely 
to prevent the onset of cooling flows. 


Meanwhile, cosmological simulations of galaxy clusters with, both, particle and grid- 
based methods have reached a fairly mature state. In particular, the addition of increasingly 
sophisticated recipes for radiative losses, star formation and stellar feedback have led to clus¬ 
ter models that can produce many observed features of galaxy clusters (Loken et al. 2002). 


Motl et al. (2004) have simulated the formation of cool cores and observed that any 
“cooling flow” is overwhelmed by the velocity held inside the cluster, which has speeds of 
up to 2000 km s“^. Nonetheless, such violent motions did not prevent the formation of cool 
cores. This shows that the formation of cool cores is inevitable unless some source of heating 
is present. Full 3D cluster collision simulations by Ritchie & Thomas (2002) and Ricker 
& Sarazin (2001) show that mergers can disrupt cooling hows. It is a common feature of 
all cosmological simulations of galaxy clusters that the IGM shows a substantial velocity 
held with many motions being supersonic. The IGM shows a complex dynamics with cool 
fronts, hlaments, shocks etc. New observations conhrm this picture. Detailed observations of 
unprecedented resolution by the latest X-ray observatories have revealed a rich portfolio of 
substructure in galaxy clusters. For example, Schuecker et al. (2001) hud substructure in the 
majority of clusters in their REFLEX-I-BGS cluster sample. The ubiquity of substructure 
points to a high frequency of mergers and other events that prevent the cluster from relaxing 
to a smooth, unperturbed state. 


Both, observations and simulations suggest that the IGM is in violent motion and that 
relaxed clusters, in which the gas sits almost statically in its potential well, are very rare, if 



4 


they exist at all. In the presence of viscosity, a fraction of the kinetic energy in these motions 
can be dissipated to heat the ICM. Thus, even in inactive phases of a central AGN, there 
can be heating in the form of viscous dissipation. 

In this paper, we study the effect of viscosity on the intracluster medium in a cosmo¬ 
logical simulation. In particular, we wish to compute to what extent viscous dissipation 
of random motions in a cluster can contribute to the heating of cluster cores. In the next 
section, we describe the technique and setup of our simulations. Finally, the results are 
discussed in Sec. 3. 


2. Simulation 

The simulations were performed with the hydrodynamics code Enzo developed by Bryan 
& Norman (1997) and Norman & Bryan (1999). Enzo is a grid-based hybrid code that cou¬ 
ples N-body computations with Eulerian hydrodynamics. It is parallelised using the MPI 
(Message Passing Interface) library, and is designed to perform simulations of cosmological 
structure formation. In Enzo, the resolution of the equations is adaptive in, both, time and 
space. In space, more hnely resolved child grids are produced where user-specihed rehnement 
criteria are met. In our simulations, each child grid has a factor 2 higher spatial resolution 
than its parent grid and a factor 8 higher mass resolution. For our purposes, the rehnement 
is controlled by the overdensity of dark matter and baryons, as well as the temperature. In 
our setting, every cell whose overdensity was higher than 4 and had a temperature above 10® 
K was hagged for rehnement. All cells on a given level are advanced with the same time step. 


Enzo uses a particle-mesh N-body method to calculate collisionless particle dynamics. 
The dark matter particles are distributed onto the grids using cloud-in-cell interpolation. 
The gravitational potential is calculated on a periodic root grid using East Fourier Trans¬ 
forms. To calculate the potential on the more hnely resolved subgrids, a multigrid relaxation 
method is used. Forces are computed on the mesh by hnite diherences, and then interpolated 
onto particle positions. 


In our simulations, we included radiative cooling, assuming that all species in the bary- 
onic gas are in equilibrium. The cooling was calculated from a tabulated cooling curve 
(Westbury & Henriksen 1992) for a plasma of hxed 0.3 solar abundance. The cooling curve 
is truncated at a minimum temperature of 10^ K. 

We used a star formation model by Cen & Ostriker (1992). According to this recipe, gas 
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is converted into collisionless star particles in regions that are Jeans unstable, contracting 
(V-n < 0) and cooling rapidly. If these three conditions are met inside a computational cell, 
a fraction of the cell’s baryonic mass is converted into a collisionless particle. 


In O’Shea et al. (2005), Enzo was compared to SPH simulations using the GADGET 
code. General agreement in the distributions of temperature, entropy and density within 
clusters was found between both numerical methods. 


Our simulations have a comoving box size of 64 Mpc/h. We chose a standard flat AGDM 
cosmology with Ob = 0.044, 0^ = 0.27, Oa = 0.73, h = 0.71 and ag = 0.9. We started from 
initial conditions at a redshift of Zinit = 60. With 8 levels of rehnement and a root grid of 
64^ cells, we achieve an effective resolution of 3.9 kpc/h. The simulations were run on 32 
processor of an IBM p690 shared-memory machine. 


The evolution of internal energy is followed by solving the energy equation 


djpe) 

dt 


+ V ■ {pe + p)v = pv ■ g pAisc, 


( 2 ) 


where e is energy per unit mass, p density, v velocity and p pressure. The dissipation of 
mechanical energy due to viscosity, per unit mass of the fluid, is given by (Batchelor 1967, 
Shu 1992, Landau & Lifshitz 1997) 


^VISC 1 

p V 3 

where A = eu and 

1 / dvi dvj\ 

2 \ dxj dxi J 

and where p is the dynamical coefficient of viscosity. In our simulations with viscosity, we 
use a third of the standard Spitzer viscosity for an unmagnetized plasma (Spitzer 1962), 
which is p = 6.0 x 10-172^5/2 

g cm ^ s taking the Goulomb logarithm to be In A = 37. 
For simplicity and want for any better model, we assumed that p is constant over cosmic 
time. This is a simplification that is unlikely to be strictly true because magnetic fields that 
inhibit viscosity will only be generated as structures evolve in the universe. It is not clear 
when and how magnetic fields form (see, e.g., Briiggen et al. (2005)). However, at early 
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times, certainly at redshifts higher than 3, viscous heating is relatively unimportant because 
the temperatures are quite small. Consequently, it appears unlikely that this simplihcation 
will affect our conclusions. 

The timestep was modihed to be the smaller value of the viscous and the Courant time 
step, i.e.: 


dt = min((itcour, dCisc) , 

where dtcour is the Courant time step and dtvisc = 0.1mm{dx^/fi). 

Velocity diffusion was simulated by solving the momentum equation 


dipvi) d , , dP di^ik 

+ j—{pvkv,) + — = + ^, 

ut uXk t)Xi uXk 


where 


'^ik 


A. 

dxi 


V 

and all other symbols have their usual meaning. 


2p, (€-ik „ 


(5) 


( 6 ) 

(7) 


3. Results and discussion 

We produced two simulations: one without viscosity and one with a third of Spitzer 
viscosity. Both runs included radiative cooling and star formation. A density projection of 
the viscous simulation is shown in Fig. 1. The largest cluster is situated near the upper edge 
of this hgure. 


Halos were identihed using the HOP-algorithm developed by Eisenstein & Hut (1997). 
The parameters of the most massive cluster in our simulation are summarised in table 1. 
The virial radius is calculated for an overdensity of dp/p = 200, and the virial mass is the 
total mass (dark matter + baryons) within the virial radius. Lx is the total X-ray lumi¬ 
nosity within i?vir in the band from 0.1 - 2.4 keV. Finally, Mstars is the total mass converted 
to stars. It is striking that, while the gross features of the cluster, such as virial radius, 
mass and temperature, are nearly identical between the runs, the X-ray luminosities differ 



substantially. The X-ray luminosity depends on the density squared and is thus dominated 
by the central portion of the clusters. As discussed in the next section, there are pronounced 
differences in the central densities, which do not affect the total mass, though. Also, the 
mass converted to stars differ between the two runs, with the viscous run producing less stars. 


In Fig. 2 - 5 we show the mass-weighted temperature, density, cooling time and entropy, 
respectively, as a function of radius in the most massive, non-merging cluster in our sam¬ 
ple. While a cool core forms in the run without viscosity, it is absent in the viscous run. 
As can be seen from the temperature prohle. Fig. 2, viscosity has essentially removed the 
cool core and the mass-weighted temperature even rises slightly in the centre. The density 
prohles also show signihcant differences (see Fig. 3). In the presence of viscosity, the den¬ 
sity is nearly constant in the core, whereas, in the non-viscous run, it rises sharply in the 
inner 40 kpc. With 1/3 Spitzer viscosity, the density is very flat over the inner hundred 
kiloparsecs. In summary, we hnd that the runs with viscous dissipation lead to a hotter and 
less dense core. Consequently, the cooling time increases in the center, with respect to the 
non-viscous runs (see Fig. 4). The central cooling time is about two orders of magnitude 
higher than in the non-viscous case and larger than the Hubble time. If viscous heating was 
this efficient, no other sources of heating would be required to prevent a cooling catastrophe. 
The entropy, which is shown in Fig. 5, displays a very extended floor and has no central 
dip as in the run without viscosity. A comparison of a statistically relevant sample of sim¬ 
ulated clusters can be compared with observed samples, such as the one by Donahue et al. 
(2005). Thus, one may be able to constrain the viscosity of the ICM from X-ray observations. 


In Fig. 6 we compare the X-ray luminosity from the cluster in the band from 0.1 - 2.4 
keV. In the case with viscosity, it is apparent that the cluster is more extended and that it 
lacks a strong emission spike in the center. 


Table 1. Properties of the most massive cluster in our sample 



without vise 

with vise 

R-virial 

1.11 Mpe 

1.12 Mpe 

Alvirial 

2.2 X 10^^ Mo 

2.3 X 10^^ Mo 

Lx 

1.6 X 10^® ergs“^ 

2.2 X 10^® ergs“^ 

Tvir 

3.0 X 10^ K 

3.1 X 10^K 

Mgtars 

7.4 X 10^2 

6.5 X 10^2 Mq 




The temperature-dependence of viscosity biases this mode of heating towards hotter, 
more massive clusters, and will, thus, affect cluster scaling relations. We hnd that the ef¬ 
fect of viscosity becomes systematically less important with decreasing mass of the cluster. 
In Fig. 7 - Fig. 9, we show the corresponding prohles of a smaller cluster with a mass of 
9.7 X 10^^ Mq. Evidently, viscosity only has a minor effect on the central density prohle. A 
proper study of this effect requires a bigger simulation box that contains a large number of 
clusters, and is the subject of future work. 


Note that in the simulations presented here, heat conduction has been neglected. As 
argued above, heat conduction and viscosity are intrinsically linked, and it would be inter¬ 
esting to study their joint effect on the cluster. Since the suppression factors of viscosity 
and conductivity are unknown, the consideration of heat conduction introduces another, es¬ 
sentially free, parameter into the problem. For this first simulation, we decided to study the 
isolated effect of viscosity only. 


In order to validate our code modules, we have repeated a low-level test simulation in 
a small box with the FLASH code. The viscosity routines in FLASH have been tested in 
Ruszkowski et ah (2004a,b). Starting from similar initial conditions and including the same 
physics, we obtain very similar results between Enzo and FLASH. 

As a result of numerical diffusion, the effective Reynolds number will be hnite, even 
in the zero-viscosity case. The effective Reynolds numbers attainable in the simulation are 
proportional to the number of grid points across the fluctuation of interest to the power n, 
where n = 3 is the order of the numerical scheme^ (Porter & Woodward 1994). As expressed 
in Eq. 1, the Reynolds number in the ICM is of the order of 50, assuming 0.3 of the Spitzer 
value. Test simulations with FLASH and Enzo suggest that the effective Reynolds number 
for the resolution chosen here is > 1000. Both runs that we have presented here have identi¬ 
cal spatial resolution. Hence, the differences between the runs are solely due to the physical 
viscosity. 


Both, radiative losses and viscosity, are sensitive to the spatial resolution of the compu¬ 
tational grid. We experimented with different rehnement criteria and compared results from 
runs with different effective resolutions. Generally, the effect of physical viscosity becomes 
larger as the resolution increases. The same is true for radiative losses. The results pre- 


^See, e.g., Bowers & Wilson (1991) for the definition of “the order of the numerical scheme”, as it is 
different from the customary definition of accuracy of a perturbative calculation. 
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sented here seem to be reasonably converged and differ only marginally from the run with a 
rehnement level less. 


Assuming hydrostatic equilibrium when inferring cluster masses from synthetic obser¬ 
vations, leads to masses that are systematically lower (by 10-15 per cent) than the actual 
cluster masses in simulations. This affects cluster constraints on cosmological parameters 
(see e.g., Allen et ah (2004)). If the ICM is viscous, it is conceivable that the assumption of 
hydrostatic equilibrium is better met. Thus, gas viscosity may reduce systematic deviations 
of cluster masses inferred from X-ray observations from their true values, which, in turn, 
may have consequences for precision measurements of cosmological parameters. 


4. Summary 

In a simple numerical experiment, we added physical viscosity to a cosmological simula¬ 
tion and studied its effect on cluster properties. The results of this experiment are interesting. 


A viscosity not too far off from the unmagnetised Spitzer value can alleviate the problem 
of cool cores in galaxy clusters. Viscosity raises the temperature in the core and reduces 
the central cooling times significantly. This effect is important for massive clusters, but 
becomes ineffective for smaller clusters and groups. For 1/3 of Spitzer viscosity, the effect 
is important for cluster above ~ 10^^ Mq. Feedback from star and AGN has been ignored 
in this work. These are additional sources of heat that will increase viscous heating. The 
interplay of various feedback processes with a possibly viscous ICM may complicate matters 
and should be looked at in future work. 


As discussed above, the actual value of the viscosity in the ICM is subject to specula¬ 
tion, and the viscosity is essentially a free parameter. Here, we assumed that the magnetic 
field is not too efficient in suppressing viscosity and took a hducial value of a third of Spitzer 
viscosity. Value cited in various studies in the recent literature assume value between 1/3 
and 1/10 of the Spitzer value. Only if magnetic helds should suppress the viscosity by 2 
orders of magnitude or more, the effect will start to be negligible. We defer the study of the 
effects of viscosity on the statistical properties of clusters and groups to a future publication 
based on a bigger simulation box containing a large number of clusters. 


Real clusters may possess complex magnetic held conhgurations that lead to a very 



inhomogeneous and anisotropic viscosity. This, in turn, can lead to complex dynamics of the 
ICM. It will be a challenge for future work to disentangle the effects of transport processes 
and magnetic helds. 
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Fig. 1.— Projected gas density at z = 0.25 showing the most massive clnster in the simnla- 
tion volnme. The side length in this hgure is 60/h Mpc. 
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Fig. 2.— Mass-weighted temperature profiles for the most massive cluster in our simulation 
(see table 1). The solid line corresponds to the run with viscosity, while the dotted line 
corresponds to the run without viscosity. 
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Fig. 3.— Density profiles for the most massive cluster. The units are Mq/Mpc“^. The 
solid line corresponds to the run with viscosity, while the dotted line corresponds to the run 
without viscosity. 
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Fig. 4.— Cooling time (in Gyrs) profiles for the most massive cluster. The solid line 
corresponds to the run with viscosity, while the dotted line corresponds to the run without 
viscosity. 
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Fig. 6.— X-ray flux of the most massive cluster normalised to the maximum. The left panel 
shows the cluster simulated with viscosity. The right panel shows the cluster simulated 
without viscosity. 
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Fig. 7.— Mass-weighted temperature profiles for a smaller cluster with a mass of 9.7 x 
10^^ Mq. The solid line corresponds to the run with viscosity, while the dotted line corre¬ 
sponds to the run without viscosity. 
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Fig. 8.— Density profiles for a smaller cluster with a mass of 9.7 x 10^^ Mq. The units 
are Mq/Mpc”^. The solid line corresponds to the run with viscosity, while the dotted line 
corresponds to the run without viscosity. 
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Fig. 9.— Cooling time (in Gyrs) profiles for a smaller cluster with a mass of 9.7 x 10^^ Mq. 
The solid line corresponds to the run with viscosity, while the dotted line corresponds to the 
run without viscosity. 




